css2<- read.csv("css05_cv.csv")
css2$part<- css2$cash + css2$vouchers
css2$part[css2$part==2]<-1
css2$part<- css2$part*100
css2$voucher_yr<-0
css2$voucher_yr[css2$elec_year>2016]<-1
css2$republican<- 0
css2$republican[css2$party=="Democratic"]<- 0
css2$republican[css2$party=="Republican"]<- 1
css2$republican[css2$party=="Non-Partisan"]<- 0
css2$race<- factor(css2$race, levels = c("European","East and South Asian", "Hispanic and Portuguese", "Likely African-American","Other"))
css2$iswhiteperson<- 0
css2$iswhiteperson[css2$race=="European"]<- 1
css2$iswhiteperson[is.na(css2$race)]<- NA
css2$income<- as.numeric(paste0(substr(css2$income,2,nchar(css2$income))))
css2$income_perc<- ceiling(ecdf(css2$income)(css2$income)*10)
css2$amount<- css2$cash_amt + css2$vouchers_amt
css3<- subset(css2, !(is.na(css2$race)))

race_ses<- data.frame(1:100)
race_ses$se<- NA
set.seed(111)
ids_to_samp<- unique(css3$id)
row_num_id<- bind_cols(css3$id,1:NROW(css3))
names(row_num_id)<-c("id","rownum")
#these take a bit of time, 
#maybe about 6? minutes to get through a 1:100 loop on my machine
#race
#estimates the coef 100 times using bootstrapped sample
#sampling at the individual level
for (i in 1:100) {
  this_it_ids<- data.frame(sample(ids_to_samp, NROW(ids_to_samp), replace = T))
  names(this_it_ids)<-"id"
  this_it_ids<- left_join(this_it_ids, row_num_id)
  css4<- css3[this_it_ids$rownum,]
  qqq1<- feols(part~iswhiteperson*voucher_yr|as.factor(id)+as.factor(elec_year) , data=css4)
  race_ses$se[i]<- qqq1$coefficients[1]
  print(i)
}

#income
css3<- subset(css2, !(is.na(css2$income_perc)))
income_ses<- data.frame(1:100)
income_ses$se<- NA
set.seed(222)
ids_to_samp<- unique(css3$id)
row_num_id<- bind_cols(css3$id,1:NROW(css3))
names(row_num_id)<-c("id","rownum")
for (i in 1:100) {
  this_it_ids<- data.frame(sample(ids_to_samp, NROW(ids_to_samp), replace = T))
  names(this_it_ids)<-"id"
  this_it_ids<- left_join(this_it_ids, row_num_id)
  css4<- css3[this_it_ids$rownum,]
  qqq1<- feols(part~income_perc*voucher_yr|as.factor(id)+as.factor(elec_year) , data=css4)
  income_ses$se[i]<- qqq1$coefficients[1]
  print(i)
}

#voting index
css3<- subset(css2, !(is.na(css2$sums)))
sums_ses<- data.frame(1:100)
sums_ses$se<- NA
set.seed(333)
ids_to_samp<- unique(css3$id)
row_num_id<- bind_cols(css3$id,1:NROW(css3))
names(row_num_id)<-c("id","rownum")
for (i in 1:100) {
  this_it_ids<- data.frame(sample(ids_to_samp, NROW(ids_to_samp), replace = T))
  names(this_it_ids)<-"id"
  this_it_ids<- left_join(this_it_ids, row_num_id)
  css4<- css3[this_it_ids$rownum,]
  qqq1<- feols(part~sums*voucher_yr|as.factor(id)+as.factor(elec_year) , data=css4)
  sums_ses$se[i]<- qqq1$coefficients[1]
  print(i)
}

#republican
css3<- subset(css2, !(is.na(css2$republican)))
republican_ses<- data.frame(1:100)
republican_ses$se<- NA
set.seed(444)
ids_to_samp<- unique(css3$id)
row_num_id<- bind_cols(css3$id,1:NROW(css3))
names(row_num_id)<-c("id","rownum")
for (i in 1:100) {
  this_it_ids<- data.frame(sample(ids_to_samp, NROW(ids_to_samp), replace = T))
  names(this_it_ids)<-"id"
  this_it_ids<- left_join(this_it_ids, row_num_id)
  css4<- css3[this_it_ids$rownum,]
  qqq1<- feols(part~republican*voucher_yr|as.factor(id)+as.factor(elec_year) , data=css4)
  republican_ses$se[i]<- qqq1$coefficients[1]
  print(i)
}

#combo reg, slightly slower than the previous four
css3<- subset(css2, !(is.na(css2$republican)))
css3<- subset(css3, !(is.na(css3$income_perc)))
css3<- subset(css3, !(is.na(css3$sums)))
css3<- subset(css3, !(is.na(css3$iswhiteperson)))
combo_ses<- data.frame(1:100)
combo_ses$income<- NA
combo_ses$republican<-NA
combo_ses$sums<-NA
combo_ses$white<-NA
set.seed(555)
ids_to_samp<- unique(css3$id)
row_num_id<- bind_cols(css3$id,1:NROW(css3))
names(row_num_id)<-c("id","rownum")
for (i in 1:100) {
  this_it_ids<- data.frame(sample(ids_to_samp, 129103, replace = T))
  names(this_it_ids)<-"id"
  this_it_ids<- left_join(this_it_ids, row_num_id)
  css4<- css3[this_it_ids$rownum,]
  qqq1<- feols(part~republican*voucher_yr+ income_perc*voucher_yr + sums*voucher_yr + iswhiteperson*voucher_yr|as.factor(id)+as.factor(elec_year) , 
               data=css4)
  combo_ses$income[i]<- qqq1$coefficients[2]
  combo_ses$republican[i]<-qqq1$coefficients[1]
  combo_ses$sums[i]<-qqq1$coefficients[3]
  combo_ses$white[i]<-qqq1$coefficients[4]
  print(i)
}

#these numbers go into the table
mean(income_ses$se) #coefficients 
mean(race_ses$se)
mean(sums_ses$se)
mean(republican_ses$se)
sqrt(var(income_ses$se)) #actual SEs
sqrt(var(race_ses$se))
sqrt(var(sums_ses$se))
sqrt(var(republican_ses$se))

mean(combo_ses$income)
sqrt(var(combo_ses$income))
mean(combo_ses$republican)
sqrt(var(combo_ses$republican))
mean(combo_ses$sums)
sqrt(var(combo_ses$sums))
mean(combo_ses$white)
sqrt(var(combo_ses$white))
